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The effect of Prandtl number on the linear stability of a plane thermal plume is analyzed under 
quasi-parallel approximation. At large Prandtl numbers {Pr > 100), we found that there is an ad- 
ditional unstable loop whose size increases with increasing Pr. The origin of this new instability 
mode is shown to be tied to the coupling of the momentum and thermal perturbation equations. 
Analyses of the perturbation kinetic energy and thermal energy suggest that the buoyancy force 
is the main source of perturbation energy at high Prandtl numbers that drives this instability. 



1. Introduction 

The classic problem of natural-convection flow above a horizontal line heat source has re- 
ceived considerable attention during the last few decades (Batchelor 1954; Fujii 1963; Gebhart 
et flZ. 1988). The temperature of the heat source is larger than that of the ambient fluid, and the 
resulting density difference creates a plume that rises up against the gravity. For steady laminar 
plumes, the similarity solutions of the pertinent boundary layer equations have been published by 
many researchers (Fujii 1963; Gebhart, Pera & Schorr 1970; Riley 1974); experimental studies 
on laminar plane plumes are in good agreement with similarity solutions (Riley 1974). Exper- 
iments (Pera & Gebhart 1971) have confirmed that the laminar plumes are unstable, and they 
sway in a plane perpendicular to the axis of the source. Pera & Gebhart (1971) have shown that 
the initial instability of plane plumes to two-dimensional disturbances can be analyzed by the 
linear stability theory and developed the coupled Orr-Sommerfeld type equations using a quasi- 
parallel flow approximation. Since Squire's theorem holds for natural convection flows (Gebhart 
et a/. 1988), it is sufficient to consider two-dimensional disturbances for the stability analysis of 
a thermal plume. 

Strictly speaking, the thermal plume is a non-parallel flow field, and the streamwise variations 
of both the laminar and disturbed flows should be incorporated in the stabiUty analysis (Hieber 
& Nash 1975; Wakitani 1985). From a weakly non-parallel spatial stability analysis (Wakitani 
1985), it has been shown that the critical Grashof number of a plane thermal plume is slightly 
larger than that predicted from the quasi-parallel theory, even though its precise value depends 
on the flow quantity (fluctuating kinetic energy or thermal energy, etc.) that is being monitored to 
calculate non-parallel corrections. It was shown that a lower branch of the neutral stability curve 
in the (frequency, Grashof number)-plane exists when the non-parallel corrections are taken into 
account. The upper branch of the neutral curve at moderate-to-large values of Grashof number 
remains relatively unaffected, however, with non-parallel corrections. 

Two non-dimensional numbers involved in natural convection phenomena are the Grashof 
number (Gr), the ratio of the buoyancy force and the viscous force, and the Prandtl number 
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(Pr), the ratio of the kinematic viscosity and the thermal diffusivity. In terms of Pr, there are two 
limiting cases: zero Prandtl number Umit (e.g. molten metals Pr ~ 10"^) and infinite Prandtl 
number limit (e.g. Pr ~ 10^ for magmas, and Pr ^ 10^^ for Earth's mantle plume). Geological 
flows involve fluids with large Prandtl numbers (Worster 1986; Lister 1987; Grossman & Lohse 
2000; Kaminski & Jaupert 2003; Majumder, Yuen & Vincent 2004) and are studied in the limit 
of infinite Prandtl number (Wang 2004) for which the inertial terms in the momentum equations 
are neglected. 

The goal of the present work is to understand the stability characteristics of high Prandtl num- 
ber plane thermal plumes. To the best of our knowledge, all stability analyses of plane thermal 
plumes (Pera & Gebhart 1971; Hieber & Nash 1975; Wakitani 1985) are confined to fliat of air 
(Pr = 0.7) and water (Pr = 6.7). We use the quasi-parallel approximation to analyse the linear 
stability of a thermal plume which is found to be unstable for very small Grashof numbers at any 
Prandtl number. At high Prandtl numbers, we find a new instabihty loop which is shown to be 
tied to the coupUng of the hydrodynamic and thermal perturbation equations. An analysis of the 
perturbation energy unveils the driving mechanism of this instabihty. 



2. Governing equations and base flow 

We consider the convective flow generated above a line heat source in an otherwise stagnant 
fluid which is maintained at a constant temperature Tc^. Let the Cartesian coordinate system is 
(x, y), with X being directed along the flow direction (i.e. against gravity) and y is the transverse 
direction, and u and v are the corresponding velocity components along x and y directions, 
respectively, and t is the time. With Boussinesq approximation, the governing equations for the 
velocity and the temperature fields are given by 

du du du „9 1 dp ^ , 

— +u—+v— = iyV\--^+gP{T-T^) (2.1) 

ot ox ay p ox 

dv dv dv „9 1 dp 

j^+uj-+v— = i^y\- -jf (2.2) 
at ox oy pay 

— + u— -h t;^ = kV^T, + = (2.3) 

dt dx dy dx dy 

Here p is the mean density of the fluid, p is the pressure, g is the acceleration due to gravity; 
the thermo-physical properties of the fluid are the thermal expansion coefficient /?, the kinematic 
viscosity v, the thermal conductivity k, the specific heat at constant pressure Cp and the thermal 
diffusivity n = k/ pcp. The boundary conditions on velocity and temperature are: u = v = 
0, T = Ts atx = y = itadu = v = T = at x'^ + y"^ ^ oo. 

2. 1 . Base flow: similarity solution 

The steady laminar base flow is given by the leading-order boundary-layer equations (Fuj ii 1 963 ; 
Gebhart et al. 1970; Pera & Gebhart 1971 ; Riley 1974) that can be expressed in terms of a stream 
function: u = dip/dy, v = —dip/dx. The resulting partial differential equations (not shown) 
can be transformed into a set of ODEs in terms of a similarity variable r] = y/5, with S = Ax/G, 
where 

Gr^iSmi^I^ „C = 4(91)"' ,2.4, 

i^^ V 4 / 

are the local Grashof number and the 'modified' Grashof number, respectively, and To(x) = 
T{x, 2/ = 0) is the local centerhne temperature. The non-dimensional stream-function and tem- 
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Figure 1 . Variations of the base-state (a) velocity (/') and (6) temperature {h) with Prandtl number, 
perature are defined via 

HV) = 77773^ — ^ /..TT i^o!-ys ^ (2-5) 



^^^^ - Z&' 



where Uc = vG^ /A:X is the local convective velocity and Tc = (To(x) — Too) = vUc/gfiS'^ is 
the local excess temperature at the centerline of the plume. With the assumption of power-law 
variation of Tc x^'^^^), the similarity equations can be obtained as 



/ + 



12 „ 4 ,2 „ 12 , ' 

.//--/ h +—Pr{fh) =0, 



(2.6) 



5 5 5 

where the prime denotes differentiation with respect to 77, and the related boundary conditions 
are (Gebhart ef fl/.1970): /(O) = /"(O) = h' (0) = 0, /i(0) = 1, /'(oo) ^ 0, h{oo) ^ 0. 
These equations have been solved by using the fourth-order Runge-Kutta method with Newton- 
Raphson correction. The far- field boundary condition was implemented at — 12, and the results 
were checked by using different values of 77 = 8, 16, 20, 50. Figure 1 shows the velocity and the 
temperature profiles for a range of Prandtl numbers. With increasing Pr, the velocity profile 
flattens across the plume width, and the temperature boundary layer becomes narrower 



3. Linear stability analysis: quasi-parallel approximation 

To analyse the stability of a thermal plume, we decompose each dynamical variable into a 
mean part (base flow) and a small-amplitude perturbation: 

u{x,y,t) ^ u{x,y) + u(x,y,t), v{x,y,t) = v{x,y) + v(x,y,t) (3.1) 
T{x,y,t) ^T{x,y) + f{x,y,t), p{x,y,t)^p{x,y)+p{x,y,t) (3.2) 

with the base flow being taken as that given by the similarity solution i2.6\ over which the per- 
turbation equations are linearized. The base flow quantities v and the x-derivatives of u and T 
are taken as zero in the linearized perturbation equations- this is called the quasi-parallel ap- 
proximation. The perturbations are assumed to be of the form such that their amplitudes depend 
on the similarity variable ry, as does the base flow. 

As in the case of the base flow, it is straightforward to show that the perturbation equations 
can be expressed in terms of a stream function: u = d^/dy, v = —dil)/dx. The resulting 
perturbation equations (not shown) are amenable to normal-mode analysis: 

(Vi,f)(a;,ry,i) = (<^,s)(r;)e'("---*), (3.3) 
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where a and u: are the non-dimensional wavenumber and frequency, respectively, with S and 
r — S/Uc being the reference length and time scales. The amplitudes for the perturbation stream 
function and temperature are made dimensionless via (j> — cj)/UcS and s — s/Tc- Substituting 
the normal-mode decomposition ( 13.31 ) into the linearized perturbation equations, the coupled 
Orr-Sommerfeld stability equations are obtained: 



2a^6' +a^d> + s ] = iaG 



s ~ o? s — iaPrG 
The boundary conditions on (j){r]) and 3(77) are: 



a/ 



(3.4) 
(3.5) 



(l){±oo) = {±00) = s(±oo) = 0. (3.6) 

3.1. Varicose and sinuous modes 

For varicose modes, both the velocity and temperatures are symmetric about the mid-plane (r] = 
0) which can be translated into following conditions on cj) and s: 

(/)(O) = 0"(O) = .s'=O, (3.7) 

whereas for sinuous modes both the velocity and temperatures are asymmetric about the mid- 
plane {rj = 0) for which the conditions on (f> and s are: 

,^'(O) = 0"'(O) = s(O) = O. (3.8) 

It is known that the thermal plumes are more unstable to sinuous modes (Pera & Gebhart 1971) 
which has been confirmed in our study too. Hence, all results are presented only for sinuous 
perturbations (I3.8I) . 

3.2. Generalized eigenvalue problem and numerical method 

The linear stability equations ( I3.4H3.5I) . along with boundary conditions ( 13.61 13.81 ). constitute a 
generalized eigenvalue problem: 

A* = AB*. (3.9) 

For the temporal stability analysis, A = a; is the eigenvalue, <i> — {</), s)'^ is the eigenfunc- 
tion and A and B are 2 x 2 matrix differential operators whose elements can be easily ob- 
tained from ( 13.4113.51 ). For the spatial stability analysis, the spatial eigenvalue X — a appears 
nonlinearly in ( 13.4113.51 ) which are subsequently transformed into a linear problem in wavenum- 
ber (a) by using the 'companion-matrix' method (Bridges & Morris 1984). For this case, <i> = 
{a^(j), 0^4), a4>, 4>, as, s)'^ is the eigenfunction and A and B are 6 x 6 matrix differential oper- 
ators; the non-zero elements of A are (with D — d/drj): An ~ —iGf, An — vjjG + 2D^, 
yli3 = iGf'D^~iGf"', Au ^ -D'^-iLoGD^, A^q = -D, A53 = iPrGh' , A55 = -iPrGf' , 
^56 = + iwPrG, A21 — 1 ~ A^2 — j443 = Ags, and B is an unit diagonal operator. 

For the temporal stability, the wavenumber, a, is real and the frequency, uj — ujr + iuJi, is 
complex, with cj, being the 'temporal' growth/decay rate of the perturbation. For the spatial 
stability, the frequency, lo, is real and the wavenumber, a — + ia^, is complex, with at 
being the 'spatial' growth/decay rate. In either case, the flow is said to be stable/unstable if uji or 
—ai >, < 0, respectively, and neutrally stable if LOi or ai ~ 0. 

For both temporal and spatial analyses, the differential eigenvalue problem ( |3.9l l is transformed 
into an matrix eigenvalue problem by discretizating the related differential operators along the 
non-periodic 77-direction. We have used two numerical methods for discretization (Malik 1990): 
(1) the finite difference method with second-order accuracy; (2) the Chebyshev spectral collo- 
cation method. The resulting matrix-eigenvalue problem has been solved by the QZ-algorithm 



a 



3 
2.5 
2 
1.5 
1 

0.5 

0, 



New instability mode in a plane thermal plume 

2 



(n) 






/ 0.005 




I 0.01 J X 

B " 




OO O-O^' ■ 






1.5 



a 



0.5 



(b) 






^^005^^^^^ 




--^01 




0.015 


\l c>0 















50 



100 150 200 250 
G 



5 
G 



10 



Figure 2. For the temporal analysis, the stabiUty diagrams in the (q, G)-plane at Pr — 200; panel b is 

the zoomed part of low-G region. 
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Figure 3. For the spatial analysis, (a) the stability diagram in the (cj, G)-plane at Pr = 200; (b) the 
variation of spatial growth-rate (— Qi) with G for = 0. 



of Matlab. The growth-rate and the phase speed of the least-stable mode, obtained from finite 
difference and spectral methods, were compared for a few test cases with different number of 
grid/collocation points (N = 101 and 151). We found that both the growth rate and the phase 
speed agreed upto the third decimal place for two methods. Moreover, from a comparison with 
published literature, we found that our neutral stability curve for air {Pr = 0.7) agrees well with 
that of Pera & Gebhart (1971). 



4. Results and discussion 

We have carried out both temporal and spatial stability analyses of a thermal plume, and most 
of the results are presented for the temporal case (except in figure[3]). 

4. 1 . Results for various Prandtl numbers 

Figure |2a) displays a typical stability diagram in the (a, G)-plane at high Prandtl numbers 
{Pr — 200) for the temporal stability analysis; figure Wib) is the zoomed part of the low-G 




region of figure|2ta). In each panel, the neutral contour (uji — 0) is marked by '0', and the flow is 
unstable inside it (see the positive growth rate contours) and stable outside. There are two distinct 
unstable zones: (a) one at low wavenumbers (a) that spans the whole range of Grashof number 
(G > Gcr), and (b) the other at relatively higher wavenumbers that spans a limited range of G. 
FigurelSJ^)) suggests that there is a minimum value of G below which the plume is stable. 

The thick line in figure |2fa-6) demarcates the regions of downstream-propagating (phase 
speed, Cr — ojr/ct > 0) and upstream-propagating (cr < 0) modes in the (a, G)-plane. The 
origin of such upstream-propagating modes (at low a) remains unclear to us at present. We have 
checked that the locus of Cr = line in the {a, G)-plane does not change by increasing the size 
of the computational domain from 77 = 12 to 77 = 100 or by increasing the number of collocation 
points. We should point out that the possibility of having upstream-propagating modes in a plane 
thermal plume (which exist for any Pr at very small values of a) has not been mentioned in 
previous works (Pera & Gebhart 1971; Hieber & Nash 1975; Wakitani 1985). Such modes might 
be analogous to certain backward-propagating modes in Ekman boundary layer (Lilly 1966) - 
this issue is relegated to a future study. 

Figure[3a) displays the analogue of figure|2a) in (uj, G)-plane for the spatial stability analy- 
sis. As expected, the stability diagram in the {uj, G)-plane also contains two unstable loops which 
are analogues of the two-loops of the temporal case, figure|2ta). Focussing on the zero-frequency 
modes (oj = 0) in figure [S^a), we plot the variation of the spatial growth rate of the least-stable 
mode (— Qfi) with G in figure [516). It is seen that the flow is unstable to w = modes beyond 
a minimum Grashof number, G ^ 0.185, for Pr — 200, and the corresponding real wavenum- 
ber is ^ 0.47. This critical point (G, a, ) = (0.185, 0.47) from the spatial analysis exactly 
matches with the intersection point between the neutral curve and the locus of Cr = modes in 
figure|2f6) for the temporal analysis. This result establishes that the upstream propagating modes 
{Cr < 0) for the temporal case are not an artifact of the numerical method. 

Here onwards, we present results only for temporal stability. Focussing on figure|2|a), we show 
the variations of the growth-rate and phase speed of the least-stable mode with wavenumber in 
figures Ufa) and|4f6) for G = 50 and G = 100, respectively. Two humps in each growth-rate 
curve correspond to two unstable loops in figure|2a), and the second hump is referred to as new 
mode since it does not have an analogue in low-Pr fluids. The discontinuities in each phase- 
speed curve correspond to crossing of different modes. 

For a range of Prandtl numbers (Pr = 0.7, 100, 200 and 500), the stability diagrams, contain- 
ing the neutral contour (oji = 0) along with a few positive growth-rate contours (iVi > 0), are 
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Figure 5. Stability diagrams at various Prandtl numbers: (a) Pr = 0.7; (b) Pr = 100; (c) Pr = 200; (d) 
Pr = 500. Thick line in each panel corresponds to the neutral contour for 'uncoupled' stability equations 
(see §4.2 for related explanation). 

compared in the (G, a)-plane in figure|3a-(i). For Pr = 0.7 (air), the stability diagram has one 
loop, and the upper branch of the neutral curve is well defined and has an inviscid asymptotic 
limit: a = 1.3847 (Pera & Gebhart 1971). In the limit G ^ oo, there exists a range of wave- 
numbers over which the flow is unstable. At high Prandtl numbers (Pr > 100), as in figure|3c), 
the neutral curve contains a kink, and there is an additional unstable loop at large a and low G. 
The size of this new unstable loop increases with increasing Prandtl number, see figure|3(i). As 
mentioned before, this new unstable loop is referred to as a new mode since it does not appear in 
low-Pr fluids. Comparing the growth rate contours for different Pr in figure|5] we find that the 
growth rate of the least-stable mode decreases with increasing Prandtl number, even though the 
size of the unstable zone in the {a, G')-plane increases in the same limit. The thick solid contour 
in each panel of figure|5]is explained in the next section. 

Now we suggest one possible experiment to realize this new instability mode. Suppose a lam- 
inar plume is disturbed by a small-amplitude sinusoidal excitation of the source with a specified 
frequency, f ~ uj/2tt (the temperature at the source is constant such that G — 100, say, in figure 
3 a). For small enough / the plume will show a wavy instability according to the lower instability 
loop in figure 3(a), however, for relatively larger / the plume is unstable to our new instability 
loop and there is a window of frequencies between two instability modes over which the plume 
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remains stable. It would be interesting to verify this transition scenario at a given Grashof num- 
ber, 'unstable^stable— >unstable' with increasing frequency, in experiments of high Pr-fluids. 

4.2. Origin of new instability loop: coupling of hydrodynamic and the nnal fluctuations 

To shed light on the origin of the new instability loop at high Pr, here we assume that the 
velocity and the temperature perturbations are decoupled from each other. The full set of stability 
equations, (|3.4H3.5l l, can be made independent from each other by dropping s and h cf) from 
equations ( 13.4b and (I3.5l l, respectively. These two sets of equations can now be solved separately 
to determine the least stable eigen-value - we have verified that the least stable mode belongs 
to the Orr-Sommerfeld equation (i.e. a purely hydrodynamic mode, eqn. I3.4l i. and the energy 
equation ( 13. 5t always yields a stable mode. 

For the uncoupled perturbations, the neutral stability curve for each Pr is superimposed as a 
thick solid contour in figure |5] (The flow is unstable inside the thick contour and stable outside.) 
A comparison of the thick line in each panel with the corresponding neutral contour of coupled 
stability equations (denoted by the thin curve 0) clearly reveals that the coupling between the 
hydrodynamic and the thermal disturbance equations is responsible for the origin of our new 
instability mode at high Pr. It is observed that the lower parts of the instability loops in fig- 
ures |3c-d) closely follow the instability loop of the 'uncoupled' Orr-Sommerfeld equation, and 
are, therefore, purely hydrodynamic in origin. 

It is clear that the coupling terms in the stability equations ( l3.4H3.5l l are responsible for ap- 
pearance of the additional instability loop at high Prandtl numbers, and solving the uncoupled 
perturbation equations would lead to incorrect results. The importance of this coupling between 
hydrodynamic and thermal perturbations at high Pr can be understood from the fact that the 
gradient of the base-flow temperature, h , (which appears in the energy perturbation equation) 
increases with Prandtl number: h'{r]) ~ Pr f {rj)h{ri) ~ Pr'^, with < e < 1, and hence cannot 
be neglected at large Pr. (From an order-of-magnitude analysis of the pertinent boundary-layer 
equations, we find e = 1/2.) 

4.3. Analysis of perturbation energy: instability mechanism 

Lastly, to understand the underlying instability mechanism, we analyse different components of 
perturbation energy. The time-evolution equations of perturbation kinetic energy and thermal en- 
ergy are obtained from ( |3.4H3.5l l by multiplying them with the corresponding complex conjugate 
quantity (/)^ and s\ respectively, and integrating them from rj — Q Xo rj = oo. Considering the 
real parts, the resulting evolution equations boil down to (Nachtsheim 1963; Gill & Davey 1969) 

AEk 
dt 



= uji I exdri = / etrKdr]+ / eyed?? + / esdrj, (4.1) 
Jo 



LUi I erdrj = / etrrdi] + / CTodri, (4.2) 
10 Jo Jo 



where 

eK 



(I P + a^l P) , etrK^af {4)r(t)i - 'Pr'i'i) , evD^G^^\4> -a^4> 

ct = I s 1^, etrT = ah (i^rSi - i^iSr) , ctd = -Pr~^G~^(^\s |^ + a^|s|^ 

and the suffixes r and i denote the real and imaginary parts, respectively. For hydrodynamic fluc- 
tuations, dSii /dt represents the rate of change of perturbation kinetic energy, EtrK = J etrK 
the rate of transfer of kinetic energy from the mean flow to perturbations via the Reynolds stress, 
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Figure 6. Distributions of (a) kinetic and (fe) thermal energies across the plume width for Pr = 0.7 (air), 
G = 100, a = 1.2923. Vertical line indicates the location of the critical layer. 
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Figure 7. Same as figure 6, but for Pr = 200, G = 100, a = 2.74645. 



EyD the rate of viscous dissipation, and Eb the rate of gain of kinetic energy through the buoy- 
ancy force. For thermal fluctuations, A£t /di is the rate of change of perturbation thermal energy, 
E'trT is the rate of gain of thermal energy from the mean temperature field and E^d is the rate 
of dissipation of thermal energy. 

The variations of different kinetic and thermal energy components across the plume width r\ 
are displayed in figures|6]and|7]for Pr = 0.7 (a = 1.2923) and 200 {a ^ 2.74645), respectively, 
at G = 100. (These two cases correspond to the neutral modes on the upper branch of the neutral 
contour in figures |5]2 and |5]:.) In each figure, the location of the corresponding critical layer 
is indicated by the vertical line. Since the instability mode is neutral (oj^ — 0), the net rate of 
gain of kinetic/thermal energy is Ex — '^i J e^rd?/ = = Et, denoted by the dotted zero- 
line in each panel. For Pr = 0.7, the kinetic energy gained by the perturbation mainly comes 
from the Reynold's stress term (EtrK) and a small amount is contributed from the perturbation 
buoyancy force (Eb)', the maximum amount of energy is dissipated at the center line (rj = 0) by 
viscous forces (Evd)- With increasing Pr, the rate of gain of kinetic energy by Reynolds stress 
becomes progressively smaller, and the buoyancy force (Eb) takes over as the main source of 
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perturbation kinetic energy (see figure 7a at Pr = 200) which balances the energy lost due to 
viscous dissipation. We conclude that at high Pr the contribution from the Reynold's stress term 
is negligible compared to the gain in kinetic energy by buoyancy force which drives instabiUty. 



5. Conclusions 

Based on a quasi-parallel stabiUty analysis of a plane thermal plume, we have uncovered a 
new instability loop at high Prandtl numbers. The origin of this new mode is shown to be tied 
to the coupling between the hydrodynamic and thermal fluctuations. The importance of this cou- 
pUng is tied to the increasing magnitude of the base-state temperature gradient with increasing 
Prandtl number. It is shown that the perturbation kinetic energy gained from the buoyancy force 
drives this instability at high Pr. The underlying instability mechanism differs from the well- 
known hydrodynamic instability mechanism for which the perturbation energy is gained from 
the mean flow via the Reynolds stress. In future, it would be interesting to analyse the effects of 
non-parallel corrections as well as the temperature-dependent transport coefficients on our new 
instabihty loop. 
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